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Influence of shower fluctuations and primary composition 
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We study the influence of shower fluctuations, and the possible presence of different nuclear species 
in the primary cosmic ray spectrum, on the experimental determination of both shower energy and 
the proton air inelastic cross section from studies of the longitudinal development of atmospheric 
showers in fluorescence experiments. We investigate the potential of track length integral and shower 
size at maximum as estimators of shower energy. We find that at very high energy (~ 10^^ — 10^*^ 
eV) the error of the total energy assignment is dominated by the dependence on the hadronic 
interaction model, and is of the order of 5%. At lower energy (~ lO^'^ — 10^* eV), the uncertainty 
of the energy determination due to the limited knowledge of the primary cosmic ray composition 
is more important. The distribution of shower is discussed as a measure of the 

proton-air cross section. Uncertainties in a possible experimental measurement of this cross section 
introduced by intrinsic shower fluctuations, the model of hadronic interactions, and the unknown 
mixture of primary nuclei in the cosmic radiation are numerically evaluated. 



PACS numbers: 96.40.Pq,96.40.-z,13.85.-t 

I. INTRODUCTION 

The fluorescence technique of ultra high energy cosmic 
ray (UHECR) detection was first explored in the pioneer- 
ing Fly's Eye detector and is currently being used in 
its successor, the high resolution Fly's Eye (HiRes) 0, 
as well as in the Pierre Auger Observatory Q that is 
currently under construction. The underlying idea is the 
detection of atmospheric nitrogen fluorescence light in- 
duced by the passage of charged particles through the 
atmosphere. The number of charged particles at depth 
X in the atmosphere, N{X), i.e. the longitudinal shower 
profile, can be extracted from data because N{X) is to a 
good approximation proportional to the amount of emit- 
ted fluorescence light. In this approximation, the total 
energy that goes into electrons and positrons (the elec- 
tromagnetic energy i?em from now on) is obtained by 
integration of the shower longitudinal profile 

Scm = acff / N{X)dX (1) 

"'0 

where aofi is the average (effective) ionization loss rate 
which is usually taken as a constant over the entire 
shower and is given by ~ 2.19 MeV/g cm"^ 4, 5]. 

The integral on the right hand side of Eq. rep- 
resents the total track length of all charged particles in 
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the shower projected onto the shower axis. Electrons and 
positrons constitute the bulk of the charged particles in a 
shower and contribute most to the production of fluores- 
cence light. In the following we neglect the contribution 
of muons and other charged particles to the production 
of fluorescence light, which is of the order of 2% (see 
discussion in 6]). 

It is generally assumed that the fluorescence rate is 
proportional to the ionization energy loss rate dE/dX, al- 
though this has been experimentally proved only to some 
extent 0, Hi • Consequently in order to estimate shower 
energy, there is in principle no need to convert the mea- 
sured fluorescence intensity first to a particle number, 
and then relate the total track length to the energy of 
the shower through Eq. The total ionization energy 
deposit can instead be obtained from the fluorescence in- 
tensity and can be used directly as an energy estimate 01 . 
However, as long as the lateral spread of shower particles 
is correctly accounted for To"], the conversion of fluores- 
cence light intensity to number of particles and then to 
energy through Eq. (Q) does not lead in principle to ob- 
servable errors mainly for two reasons: Firstly the ioniza- 
tion energy deposit depends only weakly on the particle 
energy, and secondly the shape of the energy spectrum 
of particles in an air shower changes only slowly with the 
traversed depth. Only in the very early evolution stage 
of a shower is the particle energy spectrum significantly 
harder than that at the shower maximum. The corre- 
sponding energy deposit is higher by up to a factor of 
1.5, but due to the small number of particles, the result- 
ing error in the energy estimation is negligible dl]. 

There are several additional factors, such as air pres- 
sure, density and humidity, that influence the relation of 
the fluorescence intensity to energy deposit and particle 
numbers. The discussion of these aspects, including the 
conversion of the observed light curve to a longitudinal 
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shower profile are beyond the scope of this work. 

In this article we investigate the longitudinal shower 
profile as an experiment-independent quantity and study 
its relation to the energy and mass of the primary parti- 
cle. In Sec. m the track length integral and the particle 
number at shower maximum are compared as energy es- 
timators under the assumption of an unknown cosmic ray 
composition above 10^^ eV. The model and mass depen- 
dence of the invisible energy carried by neutrinos and 
energetic muons is calculated for the QGSjet ^2 ^^'^ 
SIBYLL Ea H models of hadronic interactions. The 
mean position of the shower maximum, Xmax, and its dis- 
tribution is discussed as a measure of the primary cosmic 
ray composition and proton-air cross section in Sec. IIIII 
A summary and conclusions are given in Sec. IIVI 



II. 



SHOWER TRACK LENGTH AND ENERGY 
RECONSTRUCTION 



Using the particle number as primary observable, an 
accurate Monte Carlo calculation of the electromagnetic 
energy in a shower requires that all the contributions to 
the track length due to electrons and positrons are prop- 
erly accounted for. This implies computing the track 
length of electrons of a very wide energy range, including 
very small energies, since electrons with kinetic energy 
below 0.1 MeV typically account for ~ 10% of the electro- 
magnetic energy • A part of the energy of the primary 
particle is not converted to electromagnetic particles and 
hence not deposited as ionization energy in the atmo- 
sphere. High-energy muons deposit only a small fraction 
of their energy in the atmosphere and neutrinos escape 
detection completely. Therefore experiments have to cor- 
rect for the "unseen" energy to estimate the primary par- 
ticle energy. In a Monte Carlo simulation, once E^m is 
determined, the "unseen" energy can be calculated on a 
shower-by-shower basis as the difference to the primary 
particle energy. 

In this section we calculate the fraction of shower 
energy that is not detected by a fluorescence experi- 
ment and we study its dependence on shower energy, 
on the primary nucleus that initiates the shower, and 
on the hadronic interaction model needed to extrapolate 
to unmeasured regions of the phase space of the primary 
nucleus-air collision. We compare the resolution achieved 
with two different estimators of shower energy, namely 
the track length integral and the number of particles at 
shower maximum. For this purpose we simulate large 
samples of showers and extract their longitudinal profile. 

We use a fast hybrid simulation program that al- 
lows the simulation of longitudinal shower profiles of elec- 
trons and muons. The hybrid method consists of calcu- 
lating shower observables by a direct simulation of the 
initial part of the shower, tracking all particles of energy 
above E'thr — 0.01 E. Parameterizations of prcsimulated 
showers for all subthreshold particles are then superim- 
posed after their first interaction point is sampled. The 



sub-showers are described with parameterizations that 
give the correct average behavior, and at the same time 
describe the fluctuations in shower development of both 
electrons and muons. 

Electromagnetic cascades are simulated with a full- 
screening electromagnetic Monte Carlo at high energy 
and at low energy the longitudinal development of elec- 
tromagnetic sub-showers is calculated using Greisen's pa- 
rameterization p^ . The numerical approximation as 
given in [T^ is applied with electron-induced showers be- 
ing shifted by 0.8 radiation lengths ^^l- Greisen's for- 
mula is a good approximation to the numerical solution 
of the cascade equations with vanishing low-energy cut- 
off. 

In order to apply Eq. |^ to estimate shower energy 
from the number of particles given by our hybrid ap- 
proach, the factor Ocs in Eq. must be determined for 
our approximation of electromagnetic showers. We nor- 
malize the track length predicted by our hybrid method 
to the electromagnetic energy in photon initiated show- 
ers of energy i?om, turning photoproduction interactions 
artificially off to avoid that a fraction of Eg^ goes into a 
muonic and neutrino component. 



E, 



jN{X)dX 



(2) 



We obtain the numerical value a^ff — 2.32 McV/g cm^^, 
which we will use throughout this paper. It is important 
to realize that aoff depends on the treatment of low en- 
ergy particles due to the different kinetic energy thresh- 
olds of the Monte Carlo simulations, and therefore our re- 
sult cannot be directly compared to the numerical value 
obtained by Song et al. {ues = 2.19 MeV/g cm'^) [|| 
who performed a CORSIKA simulation with the 
threshold of 100 keV. Song et al. estimated that about 
10% of the electromagnetic energy is carried by particles 
of kinetic energy less than 100 keV, which are neither 
included in the track length simulation nor in the cal- 
culation of aeff. The value of 2.42 MeV/g cm^'^ found 
by Risse and Heck (Fig. 7 in 11]) is not in contradiction 
with our result as it refers to a simulation threshold of 
250 keV. In contrast to Q the analysis in defines acff 
as the proportionality constant between the projected 
track length of all particles above simulation threshold 
and the total calorimetric energy, including the expected 
energy deposit of the particles falling below the simula- 
tion threshold. 



A. Unseen energy in hadron-induced atmospheric 
showers 

We turn now to the study of the unseen energy in 
nucleus-induced showers. Using the hybrid method we 
have simulated showers initiated by protons, helium, car- 
bon and iron at zenith angle 9 = 45°, down to the 
approximate observation level of the HiRes and Auger 
experiments corresponding to a vertical depth of = 
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870 g/cw?. We have performed the simulations using 
two hadronic interaction models SIBYLL 2.1 [23| and 
QGSjetOl with the aim to study the influence of the 
model predictions on the amount of unseen energy in the 
shower. 

We calculate the total track length by performing the 
integral in depth of the longitudinal profile generated by 
the hybrid method. Similar to the method applied by 
air-fluorescence experiments the simulated curve is fit by 
the Gaisser-Hillas function '^Jj 

(^max-^0)/A 



Protons 



Helium 



Ngu{X)=N„ 



X-Xn 



Xr, 



Xo 



X exp 



{X — Xynax) 

A 



(3) 



where Xq is the depth of the first interaction, Xj^ax is the 
depth at which the number of electrons in the shower is 
maximal, and iVmax is the number of electrons at maxi- 
mum. First the position of the shower maximum and the 
particle number at maximum arc found by a polynomial 
fit to the shower profile near its maximum. In a subse- 
quent fit the parameters A and Xq are determined. The 
unseen energy in a shower of energy E then follows 
from 



E^j — E — Ef. 



E - Qfcff 



NGiiiX)dX . (4) 



In figure^we plot the mean unseen energy as obtained 
in showers initiated by different nuclei and compare the 
predictions of the interaction models SIBYLL 2.1 and 
QGSjetOl. The energy that is transferred to muons and 
neutrinos decreases with shower energy for all nuclei in 
both models. The main reason for this behavior is that 
as shower energy increases the average energy of charged 
pions increases as well, and in turn their probability of 
decaying into muons and neutrinos diminishes. At fixed 
energy the unseen energy is larger in showers initiated 
by heavy nuclei than in those induced by light nuclei. 
This can be understood on the basis of the superposition 
model in which a shower induced by a nucleus of A nucle- 
ons and energy E is considered as A independent proton 
showers of energy E/A, in each of which the fraction of 
unseen energy is larger than in a proton shower of energy 
E (Fig.^. In fact, the muon number in a proton shower 
scales as - E" with a = 0.86. . .0.92 O, and ap- 
plying the superposition model it can be easily seen that 
an iron shower has about 1.5 times more muons than a 
proton shower of the same energy. 

The unseen energy reaches almost a constant value at 
shower energy above 10^^ — 10^" eV, and is fairly insen- 
sitive to composition in this energy region, the relative 
difference between the unseen energy fractions among the 
different nuclei being ~ 10%. Interestingly the difference 
between the two model predictions is much bigger and 
essentially stays the same. At ultra- high energy the un- 
certainty in the missing energy assignment is dominated 
by the model dependence and not by an unknown pri- 
mary composition. 
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FIG. 1; Fraction of the total shower energy that is not seen by 
a fluorescence detector in showers initiated by proton, helium, 
carbon and iron as a function of shower energy E. 5,000 show- 
ers at zenith angle 6 = 45° were simulated for each hadronic 
interaction model: SIBYLL 2.1 (bold circles) and QGSjetOl 
(open circles). 



The unseen energy predicted by QGSjetOl is consis- 
tently larger than the corresponding energy predicted 
by SIBYLL 2.1, the reason being that the muon mul- 
tiphcity {Ef^ > 0.3 GeV) in QGSjetOl is larger than in 
SIBYLL 2.1. For example, QGSjet predicts about 35% 
more muons than SIBYLL 2.1 in proton-induced show- 
ers at 10^° eV This difference arises mainly due to 
the different multiplicities of secondaries predicted by the 
two models. The relative difference between the satura- 
tion values of the unseen energy predicted by SIBYLL 
and QGSjet is ^ 50%. This translates to a 5% uncer- 
tainty for the total energy estimate. 

Our results for the unseen energy fraction are similar 
to those obtained in Q and recently in l^^. For exam- 
ple. Song et al. find a fraction of unseen energy of about 
7% (10%) at E = 10^° eV for proton (iron) showers and 
QGSjet98 ||]. Barbosa et al. have also used SIBYLL 2.1 
in their simulations. They obtain an unseen energy frac- 
tion of about 5% and 8% for proton and iron showers at 
10^" eV simulated with the SIBYLL model J^]. There 
are differences of the order of 1-3% between these calcu- 
lations and our results which we attribute to the different 
low-energy interaction models used for the simulations, 
the approximative character of Greisen's parametrization 
for low-energy electromagnetic sub-showers, and different 
methods of calculating the track length integral. The 
latter involves extrapolating the electromagnetic shower 
component to larger atmospheric depths. 
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B. Shower energy reconstruction from the 
longitudinal shower profile 

There are several methods to reconstruct the primary 
energy of an observed air shower profile experimentally. 
Not only the total shower track length in the atmosphere 
in Eq. but also the number of particles at shower 
maximum can serve as an estimator of the shower en- 
ergy. The latter one is of particular interest for nearly 
vertical showers where only the first part of the shower 
can be observed and the uncertainty in the calculation of 
the total track length could be dominated by the extrap- 
olation of the shower profile into the unobserved region. 
In this section we study the energy resolution that is 
achieved with these two methods. 

As previously discussed, shower energy can be calcu- 
lated from the measured track length. The procedure is 
to fit a Gaisser-Hillas or other function to the observed 
longitudinal profile i22], integrate it to obtain the total 
track length, extract the energy that goes into the electro- 
magnetic component and correct it for the unseen energy 
that is estimated by Monte Carlo simulations. 

The procedure above is subject to uncertainties be- 
cause, usually, a previously calculated mean value of the 
missing energy, averaged over many simulated showers, 
is used to determine the energy of each individual event. 
Moreover, due to fluctuations in the shower longitudinal 
profile the experiments are unable to determine the type 
of primary nucleus that initiated the shower on an event- 
by-event basis. In consequence a correction for missing 
energy averaged over different primaries must be used. 
In addition the correction for unseen energy depends on 
the hadronic interaction model used to perform the sim- 
ulations. 

We estimate the energy of the shower (Et) from the 
track length obtained as indicated above in the following 
way. 



if u) Earn = {fu)a. 



off / NGiiiX)dX 




(5) 



where (/„) is the average value of the correction for 
the energy not seen. The value of (/„) depends on 
Eom and is for simplicity taken as the arithmetic mean 
over a uniform four-component mass composition, i.e. 
ifn) = (Li /«)/4 where = E^^/E are the corrections 
for unseen energy in showers initiated by different pri- 
mary nuclei. The index i corresponds to proton, helium, 
carbon and iron induced showers. For fixed primary par- 
ticle type, was obtained as the average value of Ecm/E 
in 5,000 simulated showers. 

Alternatively, shower energy can also be estimated 
from the size of the electron distribution at shower max- 
imum A^max- The relation between A^max and energy can 
be expressed as E = gNjna,x where g must be obtained 
by Monte Carlo simulations. As before, the determina- 
tion of shower energy is subject to uncertainties because 
g has to be averaged over many showers and different pri- 
maries, and is also model dependent. Then the estimated 
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FIG. 2: Average shower size at maximum normalized to 
shower energy in GeV in showers initiated by proton (bold cir- 
cles), helium (open circles), carbon (bold triangles) and iron 
(open triangles) as a function of shower energy E. 5,000 show- 
ers at zenith angle 9 = 45° were simulated for each hadronic 
interaction model: SIBYLL 2.1 (bottom panel) and QGSjetOl 
(top panel). 



energy of a shower (i?jv) follows from A'max through the 
equation. 

En = (5)Armax , (6) 

where (g) = <?i)/4 and the index i corresponds again 
to protons, helium, carbon and iron. The gi values were 
obtained as the average of E/N^^x in 5,000 simulated 
showers fixing the type of primary particle. 

Figure 13 shows the average shower size at maximum as 
a function of energy in showers initiated by different nu- 
clei. Each point represents an average over 5,000 showers. 
The results of SIBYLL 2.1 and QGSjetOl are presented. 
It is interesting to see how the dependence of A^max on 
composition changes with shower energy. In the energy 
region of 10^^ to 10^° eV the value of A^max is rather 
composition- independent, making it a good energy esti- 
mator. However, the mass dependence is of the order of 
10% in the energy range between lO" and lO^^ eV. The 
hadronic interaction model dependence for fixed energy 
and type of primary is remarkably small ~ 1 — 2%. 

For an experiment measuring showers with a steep en- 
ergy spectrum the knowledge of the event-by-event corre- 
lation, i.e. the energy resolution, is of great importance. 
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FIG. 3; Distributions of the error in shower energy determi- 
nation for two energy estimators. The results are shown for 
different nuclei and for SIBYLL 2.1. Left panels E = 10^** 
eV, right panels E = 10^'^ eV. Top panels: A'max is used as 
energy estimator. Bottom panels: the energy estimator is the 
track length integral. 



FIG. 4: Distributions of the error in shower energy determi- 
nation for two energy estimators. The results are shown for 
different nuclei and for QGSjetOl. Left panels E = 10^* eV, 
right panels E = 10^'^ eV. Top panels: A'max is used as en- 
ergy estimator in the top panels. Bottom panels: the energy 
estimator is the track length integral. 



In Fig. 131 we show the resolution in energy achieved with 
Et in Eq. ©, and En in Eq. ©, for SIBYLL 2.1 at 
two representative shower energies. We plot the relative 
difference in percentage between the estimated energy 
{Et or En) and the actual shower energy. Fig. ^ shows 
the corresponding energy resolution for QGSjetOl. The 
mean values and standard deviations of the proton and 
iron distributions shown in Figs. |3| and 0] are given in 
tables m and nil 

At energies above 10^^ eV, the energy resolution 
achieved with iVmax is similar to the one obtained when 
the track length is used as energy estimator. However the 
error in E with A^max is more asymmetric than the cor- 
responding one for the track length, due to the intrinsic 
asymmetry of the distribution of A^max at £^ = 10^^ — 10^" 
eV (see for instance Fig. 11 in reference 10). The asym- 
metry is less pronounced in the case of heavy nuclei. 

Both energy estimators produce an energy resolution 
which is not strongly dependent on composition above 
E = 10^^ eV. This can be easily understood from Figs.Q] 
and 121 in which it can be seen that at E > 10^^ eV both 
the fraction of unseen energy in the shower and A'max 
are weakly dependent on the type of primary nucleus. 
The dependence on composition is more pronounced in 
the energy range 10^^ — 10^^ eV, on average a 5 — 7% 



uncertainty in the energy determination with both meth- 
ods is introduced if the nature of the primary is unknown. 
It is interesting to notice that both methods tend to un- 
derestimate the energy for proton primaries and overesti- 
mate it for iron primaries due to the assumption of a uni- 
form four-component composition in all the energy range 
above 10^'' eV. The systematic uncertainty introduced by 
the hadronic interaction model when using A^max as esti- 
mator is about 1%, comparable to the one in the case of 
the track length integral method. 



In summary, A'max is a remarkably good energy esti- 
mator that may replace the total track length at energies 
above 10^^ eV if care is taken for the asymmetric errors 
in the distributions. An A^max-based energy determina- 
tion may be superior to the track length method when 
only a small portion of the shower around the maximum 
is detected, for instance due to the limited sensitivity or 
acceptance of the fluorescence detector. 
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TABLE I: Mean value and standard deviation of the distributions of the error corresponding to the histograms for protons in 
Figs. El and m 

Model SIBYLL 2.1 QGSjetOl 

logio(E/eV) (At) a{At) {An) a{AN) {At) a{At) {An) <T(Aiv) 



17.0 


-3.2 


6.1 


-4.4 


5.4 


-3.8 


6.8 


-4.6 


5.1 


18.0 


-1.8 


2.6 


-2.4 


3.8 


-3.0 


6.7 


-2.5 


3.8 


19.0 


-0.9 


1.7 


-0.3 


3.0 


-2.3 


2.0 


-1.6 


3.1 


20.0 


-0.2 


2.0 


-0.1 


3.2 


-0.8 


2.2 


-0.1 


3.8 



TABLE IL Mean value and standard deviation of the distributions of the error corresponding to the histograms for iron in 
Figs. 01 and H 

Model SIBYLL 2.1 QGSjetOl 

logio(E/eV) (At) a{At) {An) a(Ajv) {At) a{At) {An) a(Ajv) 



17.0 


4.1 


1.6 


4.3 


2.2 


4.1 


2.3 


4.3 


2.5 


18.0 


2.4 


0.7 


2.3 


1.7 


3.4 


2.5 


2.6 


2.2 


19.0 


1.4 


0.3 


1.1 


1.4 


3.1 


3.2 


1.0 


2.2 


20.0 


0.6 


0.2 


-1.4 


1.1 


1.1 


2.7 


-0.3 


2.0 



III. SHOWER LONGITUDINAL 
DEVELOPMENT AND P-AIR CROSS SECTION 
DETERMINATION 

The most obvious, experiment-independent observable 
characterizing the longitudinal shower profile is the depth 
of maximum, Xmax- The mean depth of maximum is a 
good measure of the composition in units of the mean log- 
arithmic mass. Indeed, within the superposition model 
one expects for a primary particle with mass number A 

{X^^) - HE /A) = {X^l) - In A, (7) 

with De being the elongation rate, a weakly energy de- 
pendent parameter. The elongation rate reflects fea- 
tures of high-energy hadron production, see for exam- 
ple 0, 01 1 making the interpretation of measurements 
model-dependent. However not only the mean depth of 
maximum carries important information about both the 
primary cosmic ray composition and the features of the 
hadronic interactions, but also its distribution. A num- 
ber of analyses using Xmax distributions to infer the pri- 
mary cosmic ray chemical composition are available in 
literature, for example H^H^H^I- In the following we 
will concentrate on the X^ax distribution of showers in 
respect to the possibility and limitations of determining 
the features of the hadronic interactions, in particular 
of measuring the high-energy proton-air inelastic cross 
section. 

The correlation of the first interaction point with the 
depth of shower maximum for proton showers was first 
employed by the Fly's Eye Collaboration in 28] . An anal- 
ysis in the context of a mixed primary composition was 



done in |25|, using the Fly's Eye data and more recently 
in "29] using HiRes data. 

The probability of having the first interaction point of 
a shower, Xi^t, at a depth greater than X is 

P(Xi„t > X) (X exp(-X/Ai„t), (8) 

with the interaction length Ajnt = (tj) /a. Here (m) is the 
mean mass of the air nuclei and a denotes the inelastic, 
particle production cross section. In case of a perfect 
correlation between Xmax and ATjnt, i.e. in case fluc- 
tuations in shower development were non-existent, one 
could use directly the exponential distribution of showers 
with large Xmax to calculate Xi^t and hence the proton- 
air cross section. However, intrinsic shower fluctuations 
modify the relation between the depth of maximum dis- 
tribution and the interaction length. This modification 
is typically expressed by a factor k with 

P(X^ax >X)=B CXp(-X/A), A = fcAint . (9) 

The factor k depends mainly on the pace of energy 
dissipation in the early stages of shower evolution. Con- 
cerning the models QGSjet and SIBYLL we have the in- 
teresting situation that, if one considers the mean depth 
of maximum, (Xmax), the mean inelasticity compensates 
to some extent the differences in the cross sections. For 
example, SIBYLL predicts the larger proton-air cross sec- 
tion but at the same time a smaller inelasticity than 
QGSjet. In contrast, the ratio A/Aint, i.e. the k factor, 
is sensitive to the inelasticity fluctuations. In general, a 
model with small fluctuations in secondary particle mul- 
tiplicity and inelasticity is characterized by a smaller k 
factor than a model with large fluctuations. Under the 
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assumption of similar fluctuations in multiplicity and in- 
elasticity, a model predicting a large average number of 
secondary particles leads to smaller overall fluctuations of 
the cumulative shower profile of the secondary particles 
and hence to a smaller k factor. Therefore the parameter 
A is very sensitive to the hadronic interaction model and 
its measurement would allow one to draw conclusions on 
general features of high-energy hadron production 30]. 
Conversely, knowing the k factor for a given model one 
can measure the proton-air cross section. 
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QGSjetOl 
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FIG. 5: Inelasticity distribution in p-air collisions ai E — 

In Fig. [SI we show the inelasticity distribution in 
proton-air interactions at lO^^eV simulated with SIBYLL 
2.1 and QGSjetOl. Here inelasticity K is defined as 
K = {E — E\cad)/E with i^iead being the energy of the 
fastest baryon (p, n, or A) and E the projectile energy. 
Although the mean inelasticities predicted by QGSjet 
and SIBYLL are similar (0.77 and 0.72 respectively), 
the pronounced peaks at small and large inelasticities in 
QGSjet events induce somewhat larger shower-to-shower 
fiuctuations (see below). The large difference in the pre- 
dicted mean charged particle multiplicities (QGSjet: 540, 
SIBYLL: 315 a.t E ^ eV [H) is of lesser importance 
since most of the secondary particles are very slow. 

Fig.l^shows the Xmax distributions for proton showers 
at fixed energies as predicted by QGSjetOl and SIBYLL 
2.1. The different slopes of the tails of the distributions 
stem from the different slopes of the exponential first in- 
teraction probability (Eq. ||SJ)), and the different intrinsic 
shower fluctuations predicted by the models. In table IIIII 
we give the numerical values of the slope of the Xmax 
distribution. A, obtained by doing an exponential fit to 
the tail of the SIBYLL and QGSjet distributions using 
Eq. ©. 

In the absence of internal fluctuations, all showers 
would develop through the same amount of matter, AX, 
between the first interaction point and maximum. As 
a consequence, a perfect correlation between X^ax and 
A^int would exist, and their distributions would have 
exactly the same shape, shifted by a constant AX = 
A'max — A^int- In that case the slope of the A"max distri- 



bution A would be equal to the mean interaction length 
Aint- Table mil compares the predictions of SIBYLL 2.1 
and QGSjetOl on these three quantities. Also shown in 
the table is the standard deviation of the distribution in 
AX, which gives an idea of the size of the fiuctuations in 
the shower longitudinal profile. 

The effect of fluctuations in AX is to broaden the 
correlation of X^ax with Xint and to change its slope. 
SIBYLL 2.1 predicts less fluctuations than QGSjetOl, 
the difference between the widths of the AX distribu- 
tions is however fairly small ^ 6 — 7 g/cm^. The dif- 
ferent fluctuations of the two models are also reflected in 
the larger width of the QGSjetOl Xmax distribution com- 
pared to SIBYLL 2.1. QGSjetOl also predicts a smaller 
p-air cross section (larger Ai„t) than SIBYLL 2.1. This, 
added to the fact that the intrinsic fluctuations in shower 
development are larger in QGSjetOl, makes the slope of 
the QGSjet distribution in Xmax flatter than the SIBYLL 
one, as can be seen in Fig. El and in table IIIII Interest- 
ingly, the k factors in SIBYLL and QGSjetOl are very 
similar within their statistical errors, somewhat larger in 
QGSjetOl (table mj. This means that the difference in 
the slopes of the Xmax distributions is dominated by the 
different p-air interaction lengths predicted by the mod- 
els, implying that the larger intrinsic shower fluctuations 
of QGSjetOl play a less important role. This conclu- 
sion is different when the same analysis is done for the 
old version of QGSjet, namely QGSjet98. QGSjet98 pre- 
dicts larger fluctuations than QGSjetOl, mainly due to 
a larger diffractive cross section, added to the fact that 
both versions have the same total and inelastic cross sec- 
tion and hence the same p-air interaction length. The 
larger multiplicity predicted by QGSjetOl further reduces 
the fluctuations. As a consequence QGSjet98 predicts 
larger fc-factors than QGSjetOl and SIBYLL 2.1. Numer- 
ical values of k in QGSjet98 for energies E = 10^^, 10^^ 
and 10^0 eV are k = 1.20 ± 0.02, 1.24 ± 0.02 and 
1.16 ± 0.02 respectively, the corresponding reduced 
are xVdof = 0.56, 0.91 and 1.18 31]. 

Finally, within a single model (QGSjetOl or SIBYLL 
2.1), the fc-factors depend very weakly on primary energy. 
This is just reflecting the weak energy dependence of the 
intrinsic shower fluctuations as can be demonstrated by 
looking at the values of i7(AX) in table lllll 



Influence of fitting range on fc-factor 
determination 



There are a number of complications making the mea- 
surement of the parameter A difficult. First, the X^ax 
distribution is not a perfect single exponential. As a 
consequence the fc-factor depends on which part of the 
ATmax distribution is used for fltting. This is illustrated 
in Fig. 13 where a non-constant k factor as a function 
of the smallest Xmax considered in the flt (X^^'^), is 
shown for SIBYLL and QGSjet respectively. To test 
how well a single exponential would describe the tail of 
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TABLE III: p-air interaction length (Aint), slope of the fitted tail of the Xmax distribution (A), and standard deviation of the 
A'niax and AX = Xmax — Xint distributions (all in g/cm^), where Xmax is the depth of shower maximum, and Xint is the depth 
at which the first p-air interaction occurs. These quantities are shown for different shower energies as predicted by SIBYLL 
2.1 and QGSjetOl. 30,000 proton showers were simulated to make each of the distributions. The fit to the tail of the Xmax 
distribution in order to obtain the numerical values of A and k — A/ Aint was performed using only the trailing edge of the 
distribution, 100 g/cm^ beyond the peak of the distribution i.e., Xmax > X^^^ + 100 g/cm^. Also shown is the x^/dof of the 
fit. Errors, where shown, are statistical. 
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E [eV] 


Aint 


A 


O" ( A^max ) 


a{AX) 


k 


xVdof 


Aint 


A 


0"(Aniax) 
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k 
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10^« 


43.64 


50.63 ± 0.70 


59.03 
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58.33 ±0.87 
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39.49 


47.12 ±0.68 
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1.16 ±0.02 


1.08 


44.93 


53.58 ±0.85 


63.32 
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1.18 ±0.02 
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10^0 


35.93 


42.49 ± 0.73 
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1.14 ±0.03 


0.94 
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FIG. 6: Distribution of the depth of maximum of proton- 
induced showers at energies E = lO'^^, 10^® and 10^" eV 
(from top to bottom panels). 30,000 showers were simulated 
with QGSjetOl (solid histogram) and SIBYLL2.1 (dotted his- 
togram) . 



the Xinax distributions in both models, we made fits of 
the k factor as function of X^^^^, assuming k being con- 
stant. The mean values are k = 1.15 ± 0.01, 1.16 ± 0.01 
and 1.16 ± 0.01 at £: = lO^^, lO^^ and lO^" eV re- 



spectively for SIBYLL 2.1, the corresponding values 
per degree of freedom being x^/dof = 0.25, 0.71 and 
1.6. For QGSjetOl the corresponding values of k and 
xVdof are: A; = 1.18 ± 0.01, 1.18±0.01, 1.16±0.01and 
= 0.65, 0.60, 1.48. 

The hypothesis of a flat behavior of k versus X^^^ is 
not as bad as Fig. [7| might indicate. The large errors 
of the points at large X^^^, although they deviate most 
from a constant value, have the smallest weights in the 
fit and do not affect much x^/dof. The large errors stem 
from the lack of statistics in the far tail of the X,nax dis- 
tribution [33 |. These errors might introduce large uncer- 
tainties in the determination of the p-air inelastic cross 
section, especially if a cut at large Xmax has to be ap- 
plied in order to avoid contamination by other nuclear 
species that might be present in the primary cosmic ray 
spectrum. This last issue is the subject of subsection C 
below. 



B. Influence of resolution in Amax on fc-factor 
determination 

So far we have assumed an ideal experiment which is 
able to measure the depth of maximum with infinite res- 
olution i.e., AA"iiiax = 0. However in the real world the 
accuracy is not infinite, in fact the HiRes collaboration 
has published a value of AXmax ~ 35 g/cm^ for the res- 
olution in the depth of maximum |27l l29l|. The purpose 
of this subsection is to explore the effect of the finite res- 
olution on the numerical value of the k factor. For this 
purpose we take the A"max distribution obtained before 
assuming a perfect resolution, and we smear it with a 
Gaussian of standard deviation equal to the ATmax reso- 
lution reported by HiRes. An example of the effect of this 
smearing on the original distribution is shown in Fig. |S1 
As expected the smeared distribution is wider than the 
original one (by about 15 g/cm^). However the slope of 
the distribution is not significantly modified, and as a 
consequence the change in the k factor is still within its 
estimated statistical error. This is demonstrated in table 
IIVI where the k factors obtained from the non-smeared 



9 



E=10'°eV 



1.5 
1.45 

1.4 
1.35 

1.3 
1.25 

1.2 
1.15 

1.1 
1.05 
1 



QGSjetOl A 
SIBYLL2.1 * 



1.5 
1.45 

1.4 
1.35 

1.3 
1.25 

1.2 
1.15 

1.1 
1.05 
1 




700 750 800 850 900 950 1000 1050 1100 



QGSjetOl O 
SIBYLL2.1 • 



"T — ' — I — ' — r 



700 750 800 850 900 950 1000 1050 1100 



10000 



1000 




1400 



FIG. 8: Distribution of the depth of maximum of proton- 
induced showers at enerey E = 10^** eV. 30,000 showers were 
simulated with SIBYLL2.1. The solid histogram is the non- 
smeared, perfect Xmax resolution distribution, and the dashed 
histogram is the solid distribution after smearing it with a 
Gaussian of standard deviation AXmax = 35 g/cm^. 
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FIG. 7: Numerical values of the k factor for proton- 
initiated showers simulated with SIBYLL 2.1 (full symbols) 
and QGSjetOl (empty symbols) in dependence on the consid- 
ered minimal atmospheric depth, X^^^, above which the fit 
to the tail of the distribution in Xmax is performed. 30,000 
showers were simulated to make the distributions. 



and the smeared distributions are presented for SIBYLL 
and QGSjet. 



Influence of composition on fc-factor 
determination 



A second complication in the determination of the in- 
elastic p-air cross section arises due to the mass com- 
position of the primary cosmic rays. Contamination of 
the proton spectrum by heavier elements may lead to 
changes in the measured parameter A and hence a mis- 
interpretation of the data in terms of the cross section 
or k factor. To investigate this issue we simulated iron 
showers and contaminated the proton spectrum assum- 
ing the primary cosmic ray composition reported by the 
HiRes collaboration consisting on 70% protons and 30% 



iron [231 ■ We assumed that the composition is energy- 
independent in the energy range between 10^^ eV and 
10^° eV, as indicated by the HiRes data [27l |. 

Figs. El El show the distributions of Amax for the 
measured HiRes composition, as predicted by the mod- 
els SIBYLL 2.1 and QGSjetOl respectively. A Gaus- 
sian Ainax resolution of standard deviation AAmax = 
35 g/cm^ was folded in both the proton and the iron dis- 
tributions. The figures also show the individual contribu- 
tions. As is clearly seen, a fraction of 30% of iron nuclei 
does not significantly contribute to the total distribution 
beyond the peak. In fact the change of slope with respect 
to a pure protonic composition is very small, almost neg- 
ligible in the tail of the distribution 100 g/cm^ beyond 
the peaks of the distributions. This conclusion applies for 
both SIBYLL and QGSjet. As a consequence the k fac- 
tors do not change with respect to those given in t able II V I 
for protons with a resolution in Amax of 35 g/cm^. 

Clearly if helium, being the nuclear species that pro- 
duces the largest average Amax, is present in the cos- 
mic ray spectrum in this energy range, the contamina- 
tion might be more important and a cut well beyond the 
peak of the distribution, in the far tail of the total dis- 
tribution, has to be performed in order to avoid a bias 
in the cross section determination. The exact position of 
the cut depends on the fraction of helium in the primary 
cosmic ray spectrum. To further investigate this point we 
also simulated helium-induced showers at E ^ 10^* eV 
and we plot the distribution in A^ax for a composition 
with varying fractions of protons and helium and a con- 
stant fraction of 30% iron. This is shown in Fig.^] The 
/c-factors obtained when fitting the three Amax distribu- 
tions in fig. II II using only the region 100 g/cm^ beyond 
its peak, are: fc=1.19±0.02, 1.15 ± 0.02 and 1.07± 0.02 
for fractions of proton+heUum 70% + 0%, 50% + 20% 



10 



TABLE IV: Standard deviation of the AX = Xmax — Xint distribution, k factors and reduced of the fits performed in order 
to obtain k. These quantities are shown for a non-smeared perfect Xmax resolution distribution, and for the same distribution 
after smearing it with a Gaussian error AXmax = 35 g/cm^. The results of SIBYLL 2.1 and QGSjetOl hadronic interaction 
models are shown. Errors, where shown, are statistical. Only the trailing edge of the distribution 100 g/cm^ beyond shower 
maximum is used to make the fits. 



Model 


SIBYLL 2.1 


QGSjetOl 




g/cm^ 


35 g/cm^ 


g/cm^ 


35 g/cm^ 


E [eV] 


k xVdof 


k xVdof 


k xVdof 


k xVdof 


10i« 
10^'' 
10^° 


1.15 ±0.03 1.24 

1.16 ±0.02 1.08 
1.14 ±0.03 0.94 


1.18 ±0.03 1.80 
1.18 ±0.02 1.06 
1.20 ±0.03 1.62 


1.18 ±0.03 1.83 
1.18 ±0.02 0.90 
1.14 ±0.02 0.79 


1.19 ±0.02 0.76 
1.17 ±0.02 0.75 
1.15 ±0.02 0.50 




FIG. 9: Depth of maximum distributions for a composition FIG. 10: Same as Fig. Elfor QGSjetOl. 

consistent on 70% protons and 30% iron nuclei. From top to 
bottom panels the primary energy is _E = 10^*, 10^^ and 10^'' 
eV. The distributions were obtained using the SIBYLL 2.1 

hadronic generator. A X^ax Gaussian resolution of 35 g/cm^ and 30% + 40% respectively, keeping the fraction of iron 
was folded in. constant. Clearly, if no cut is applied at X^ax larger than 

the nominal value of 100 g/cm^ beyond the peak of the 
distribution, an important systematic bias is introduced. 
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FIG. 11: Depth of maximum distributions for a composition 
consisting of 70% protons and 30% iron nuclei (top panel), 
50% protons, 20% helium and 30% iron nuclei (middle panel), 
and 30% protons, 40% helium and 30% iron nuclei (bottom 
panel). The primary energy is E = 10^*. The distributions 
were obtained using the SIBYLL 2.1 hadronic generator. A 
Xmax Gaussian resolution of 35 g/cm^ was folded in. 

This will induce an error in the determination of the cross 
section if there is a relatively large (> 20%) fraction of 
helium nuclei in the primary cosmic ray spectrum. 

IV. SUMMARY AND CONCLUSIONS 

We have used a hybrid simulation scheme jl^ for a 
quantitative evaluation of certain systematic uncertain- 
ties in the interpretation of fluorescence measurements 
of giant air showers. The sources of uncertainty we in- 
vestigate include the model of hadronic interactions used 
for shower simulation, the unknown mixture of primary 
nuclei in the cosmic radiation and intrinsic fluctuations 
in shower development. We do not investigate uncer- 



tainties in shower reconstruction, detector acceptance or 
other technical aspects related to properties of the envi- 
ronment or performance of the detectors. 

As an illustration of uncertainties arising from the need 
to extrapolate hadronic interaction models outside the 
kinematical region and energies explored by accelerator 
experiments, we compared two specific interaction mod- 
els, QGSjet [13 and SIBYLL both of which agree 
with each other and with a range of accelerator data for 
i/s ~ TeV and below. We find that the correction for un- 
seen energy (i.e. energy lost to neutrinos and muons that 
reach the ground) is consistently larger for QGSjet than 
for SIBYLL. The difference is such that, for a given track 
length integral, the assigned energy will be about 5% 
higher when the same data are interpreted with QGSjet 
rather than with SIBYLL. 

For primary energies below about 10^^ eV/nucleus the 
fraction of unseen energy depends significantly on the 
mass of the primary nucleus. For example, at 10^^ eV 
with SIBYLL about 7% of the primary energy is lost to 
neutrinos and muons as compared to 13% for iron. 

In view of the steep cosmic-ray energy spectrum, 
knowledge of the energy resolution is of great impor- 
tance. The track length integral can be used to assign 
energy when a sufficient portion of the profile is mea- 
sured to fix the parameters needed to complete the in- 
tegral. It has an intrinsic resolution of 2-4%, depending 
somewhat on interaction model and energy (narrower at 
higher energy). Size of shower at maximum gives only 
a marginally broader energy resolution (3-5%) and can 
be used when much of the profile after maximum is not 
measured, provided care is taken to correct for a slightly 
asymmetric distribution. There is in both methods some 
dependence on primary mass of the relation between the 
measured quantity and the primary energy which leads 
to a ~ 5% systematic uncertainty if the primary mass is 
not separately determined. 

Intrinsic fluctuations in shower development (after the 
first interaction) affect the relation between the inter- 
action length (Aint) and the slope A that describes the 
exponential tail of the A'^ax distribution. The relation 
is often expressed with a 'fc factor as A = fc x Aint. Dif- 
ferences in k factors for the range of models studied here 
are at the level of 5-7%, implying a similar uncertainty 
in the p-air cross section that may be inferred from mea- 
surements of shower profiles. Further uncertainties arise 
to the extent that an unknown fraction of helium and 
other nuclei contaminate the tail of the measured Amax 
distribution. 
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